library(texreg)

#Function For Cluster Robust SEs
fit.w.robust <- function(model,cluster.var){
	robust.se <- function(model, cluster){
		require(sandwich)
		require(lmtest)
		M <- length(unique(cluster))
		N <- length(cluster)
		K <- model$rank
		dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
		uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
		rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
		rcse.se <- coeftest(model, rcse.cov)
		return(list(rcse.cov, rcse.se))
	}
	out <- model$model
	clustervar <- mapply(paste,cluster.var,out[,grep(pattern=cluster.var, x=names(out))],sep="")
	vcov <- robust.se(model, clustervar)
	model$se <- vcov[[1]]
	model$coeftest <- vcov[[2]]
	return(model)
}

load('nes.placebo.RData')

##############################################
#############
#Within-District
within.district <- lm(I(congress.vote=='Dem') ~ cov2c + I(PID=='Dem') + cov2c*I(PID=='Dem') + factor(district.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income) + lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema, data=nes.placebo[which(nes.placebo$runopposed==0),])
within.district <- fit.w.robust(model=within.district,cluster.var="district.year")

#############
#Redistricting
redistricting <- lm(I(congress.vote=='Dem') ~ cov2c_92 + I(PID=='Dem') + cov2c_92*I(PID=='Dem') + factor(fips) + factor(state.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income) + I(incumbent.party=="Dem") + age_rep + freshman + tenure + died + retired + powercom + status + incrun + I(incumbent.party=="Dem")*age_rep + I(incumbent.party=="Dem")*freshman + I(incumbent.party=="Dem")*tenure + I(incumbent.party=="Dem")*died + I(incumbent.party=="Dem")*retired + I(incumbent.party=="Dem")*powercom + I(incumbent.party=="Dem")*status + I(incumbent.party=="Dem")*incrun, data=nes.placebo[which(nes.placebo$runopposed==0),])
redistricting <- fit.w.robust(model=redistricting,cluster.var="fips")

##############################################
#PLACEBO - SENATE VOTE
within.district.senate <- lm(I(senate.vote=='Dem') ~ cov2c + I(PID=='Dem') + cov2c*I(PID=='Dem') + factor(district.year) + factor(district.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income) + factor(race) + lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema, data=nes.placebo[which(nes.placebo$senate.race.competitive==1),])
within.district.senate <- fit.w.robust(model=within.district.senate,cluster.var="district.year")

redistricting.senate <- lm(I(senate.vote=='Dem') ~ cov2c_92 + I(PID=='Dem') + cov2c_92*I(PID=='Dem') + factor(fips) + factor(state.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income), data=nes.placebo[which(nes.placebo$senate.race.competitive==1),])
redistricting.senate <- fit.w.robust(model=redistricting.senate,cluster.var="fips")

##############################################
#PLACEBO - PRESIDENT VOTE
within.district.president <- lm(I(president.vote=='Dem') ~ cov2c + I(PID=='Dem') + cov2c*I(PID=='Dem') + factor(district.year) + factor(district.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income) + factor(race) + lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema, data=nes.placebo)
within.district.president <- fit.w.robust(model=within.district.president,cluster.var="district.year")

redistricting.president <- lm(I(president.vote=='Dem') ~ cov2c_92 + I(PID=='Dem') + cov2c_92*I(PID=='Dem') + factor(fips) + factor(state.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income), data=nes.placebo)
redistricting.president <- fit.w.robust(model=redistricting.president,cluster.var="fips")

##############################################
#PLACEBO - GOVERNOR VOTE
within.district.governor <- lm(I(governor.vote=='Dem') ~ cov2c + I(PID=='Dem') + cov2c*I(PID=='Dem') + factor(district.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income) + factor(race) + lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema, data=nes.placebo)
within.district.governor <- fit.w.robust(model=within.district.governor,cluster.var="district.year")

redistricting.governor <- lm(I(governor.vote=='Dem') ~ cov2c_92 + I(PID=='Dem') + cov2c_92*I(PID=='Dem') + factor(fips) + factor(state.year) + factor(age) + factor(gender) + factor(education) + factor(urbanism) + factor(income), data=nes.placebo)
redistricting.governor <- fit.w.robust(model=redistricting.governor,cluster.var="fips")

#############################################################################################
#Output Results 

#Within District Specification
within.district.se <- sqrt(diag(within.district$se))
within.district.senate.se <- sqrt(diag(within.district.senate$se))
within.district.governor.se <- sqrt(diag(within.district.governor$se))
within.district.president.se <- sqrt(diag(within.district.president$se))

within.district.pvalue <- coeftest(within.district, vcov=within.district$se)[,4]
within.district.senate.pvalue <- coeftest(within.district.senate, vcov=within.district.senate$se)[,4]
within.district.governor.pvalue <- coeftest(within.district.governor, vcov=within.district.governor$se)[,4]
within.district.president.pvalue <- coeftest(within.district.president, vcov=within.district.president$se)[,4]

#Redistrict Specification
redistricting.se <- sqrt(diag(redistricting$se))
redistricting.senate.se <- sqrt(diag(redistricting.senate$se))
redistricting.governor.se <- sqrt(diag(redistricting.governor$se))
redistricting.president.se <- sqrt(diag(redistricting.president$se))

redistricting.pvalue <- coeftest(redistricting, vcov=redistricting$se)[,4]
redistricting.senate.pvalue <- coeftest(redistricting.senate, vcov=redistricting.senate$se)[,4]
redistricting.governor.pvalue <- coeftest(redistricting.governor, vcov=redistricting.governor$se)[,4]
redistricting.president.pvalue <- coeftest(redistricting.president, vcov=redistricting.president$se)[,4]

#Within-District
texreg(list(within.district, within.district.senate,within.district.governor,within.district.president), omit.coef=c('factor|lpop|ldensity|ld|urban|curb|school|linc|s20l|s65m|black|fema|(Intercept)'), custom.model.names=c("House","Senator","Governor","President"), stars=c(0.05), caption="Pr(Vote for Politician) by Level of Congruence - Within District Specification", override.se=list(within.district.se, within.district.senate.se, within.district.governor.se, within.district.president.se), override.pval=list(within.district.pvalue, within.district.senate.pvalue, within.district.governor.pvalue, within.district.president.pvalue))

#Redistricting Specification
texreg(list(redistricting, redistricting.senate, redistricting.governor, redistricting.president), custom.model.names=c("House","Senator","Governor","President"), stars=c(0.05), caption="Pr(Vote for Politician) by Level of Congruence - Redistricting Specification", omit.coef=c('factor|runopposed|incrun|close|chair|rankmem|leader|status|powercom|retired|died|majparty|tenure|freshman|rep|(Intercept)|incumbent.party'),override.se=list(redistricting.se, redistricting.senate.se, redistricting.governor.se, redistricting.president.se), override.pval=list(redistricting.pvalue, redistricting.senate.pvalue, redistricting.governor.pvalue, redistricting.president.pvalue))